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A straight forward numerical technique, based on the Gross-Pitaevskii equation, is used to 
generate a self-consistent description of thermally-excited states of a dilute boson gas. The process 
of evaporative cooling is then modelled by following the time evolution of the system using the same 
equation. It is shown that the subsequent rethermalisation of the thermally-excited state produces 
a cooler coherent condensate. Other results presented show that trapping vortex states with the 
ground state may be possible in a two-dimensional experimental environment. 
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I. INTRODUCTION 



Demonstrations of Bose-Einstein condensation (BEC) 
in trapped, dilute alkali gases [|l]-[| provide a strong in- 
centive to develop accurate theoretical models of the phe- 
nomenon. In particular, theories that go beyond the sim- 
plest approach based on the Gross-Pitaevskii (GP) equa- 
tion [Q are under intense study Ij-gJ. In this paper, we 
show how the GP equation can handle the time evolu- 
tion of Bose condensed gases at Unite temperature or in 
non-equilibrium states. 

Since current experiments produce BECs at a fi- 
nite temperature, the resulting Bose-condensed vapour 
is, therefore, a mixture of condensate and excitations. 
Understanding the excited states through the Gross- 
Pitaevskii equation however, has generally been lim- 
ited to studying the linearised excitations of the sys- 
tem |9|-|l^]. It is possible, using hnite temperature field 
theory, to calculate the energy and mean field due to 
the excited quasiparticles in a self consistent manner for 
static configurations |p|. To study the time evolution 
of thermally-excited condensates using such techniques 
does however, require a prohibitive amount of compu- 
tational time even for a one-dimensional system. The 
result is that there is no computationally simple way to 
study the evolution of thermally-excited states or self- 
consistent combinations of them that employs the fuller 
quantum field theory approach. This is, of course, a very 
great limitation in modelling experiments. 

The evaporative cooling process |?]]|], where hotter 
atoms are removed and the resulting rethermalisation of 
the atomic gas produces an overall cooling is, of course, 
fundamental to producing a BEC. Current theoretical 
descriptions 14 - P| generally rely on the use of kinetic 
equations or approximations to them. To date, evapora- 
tive cooling has not been studied using the GP equation, 
although it is implicit in the work of Kagan and Svis- 
tunov PQ-pa]. We shall not be studying the problems 



associated with the critical point and we refer the reader 
to the recent work of Stoof, Ref . |^3| , for a discussion of 
related issues. 

This paper addresses the evolution of a thermally- 
excited boson gas, including evaporative cooling, through 
a Gross-Pitaevskii formalism. Section |n] contains an out- 
line of the theoretical and numerical description of a 
Bose-condensed atomic vapour based on the GP equa- 
tion which is used throughout the following sections. 

In section III a numerically-efficient method for gener- 
ating self-consistent thermally-excited states of a boson 
gas in three spatial dimensions is presented. It is based on 
the GP equation which describes the thermalisation pro- 
cess or self-interaction of the mean field through the non- 
linearity. Results of a two-dimensional numerical simula- 
tion using parameters corresponding to a sample of 87 Rb 
atoms indicate that the method is robust and largely in- 
dependent of the initial distribution 

Finally, in section [V , a model of the evaporative cool- 
ing process based on the mean field description is pre- 
sented. In the model the most energetic atoms are re- 
moved from the atomic vapour, where neither the as- 
sumption of sufficient ergodicity nor approximations re- 
lating to the strength of the interactions between the 
thermal and condensate parts of the thermally-excited 
state are necessary. Results from a numerical simulation 
of this model using a similar sample of atoms are also 
included. 



II. THEORETICAL AND NUMERICAL 
BACKGROUND 



The evolution of the mean field, ^(r, t), is described 
by the Gross-Pitaevskii equation || 
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where m is the mass of the atoms and N is the number 
of atoms in the atomic vapour. The mean field is the en- 
semble average of the boson field operator, ^(r, t), where 
the nature of the ensemble will be addressed below. The 
harmonic trapping potential and nonlinear constant are 
defined by 



V(r) = rnuitr ■ r/2 
Uq = 4ir?i 2 a/m 



(2a) 
(2b) 



where uj t is the angular frequency of the trapping poten- 
tial and a is the inter-atomic s-wave scattering length. 
The self-interaction of the mean field, describing the two 
particle scattering of the atoms, is determined by the 
nonlinearity in Eq. (Q) and is characterised by the num- 
ber of atoms and their relative scattering strengths. The 
explicit appearance of the number of atoms in Eq. (Q) is 
associated with the normalisation condition 



|*(r,t)| 2 dr = 1 



(3) 



To facilitate a computational solution, Eq. (|l]) is scaled 
in terms of harmonic oscillator units 
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where C = (x, Hi z ) and e is ratio of trapping potentials. 

In a 'Bose-Einstein membrane' the atomic de Broglie 
wave is confined strongly in one-dimension, so the ratio 
of the trapping potential, e, is large. For this two di- 
mensional, large e, scenario there will be no excitations 
in the third, confined, dimension due to the much larger 
energy that they would require. Furthermore, in this con- 
fined dimension the number of atoms per unit length will 
be small allowing the mean field to be represented by a 
ground simple harmonic oscillator (SHO) state, namely 
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#(C,t) = 9(x,y,r) exp j 
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Substituting Eq. (o) into the GP equation results in 



+r/(z)|*(|,r)| 2 exp 



(7) 



c i/a-2 



where £ = (x,y). Here the nonlinear term, r/(z) = 
fi~jfejT72 1 is parameterised by the confined dimension. 



Choosing z = results in the following two dimensional 
GP equation, 

id T V(tr)= (-V 2 + ^+r?|*(£,T)| 2 W^), (8) 



where rj = 77(0) = h ^, U ^y/2 ■ Here ^(^jt) is the mean 
field surface density, to reconstruct the full volume den- 
sity one simply multiplies the surface density by the 
ground SHO state. 

The formal solution of Eq. (0) is given by 
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(9) 



(-v 2 + e 2 /4 



where St is the time step and H 
V\*(tr)\ 2 ). 

Eq. ^ allows a numerical solution using the split-step 
operator method. The solution, performed in two dimen- 
sions, exhibits features that are extremely difficult to see 
in the one-dimensional models; it also enables the study 
of the crucial process of vortex formation. Furthermore, 
this solution is a direct simulation of a two-dimensional 
'Bose-Einstein membrane', realisation of which in two- 
dimensional atomic waveguides is a subject of increas- 
ing interest, with recent experimental and theoretical ad- 
vances 

It is pertinent to note that the analysis leading to 
Eq. (^|) represents an approximation of the available 
quasiparticle modes. This approximation results in the 
phenomena presented occurring over longer time scales 
than would be expected in an experimental configuration. 
Furthermore, the quantitative time scales for the two- 
dimensional and three-dimensional cases are, of course, 
quite distinct as the corresponding increase in mode den- 
sity results in shorter time scales over which the phenom- 
ena presented occur. Nevertheless the topics presented in 
this paper (self-consistent thermally-excited state gener- 
ation, evaporative cooling and the mechanisms relating 
to them) are not specifically related to the number of di- 
mensions considered. Therefore, the arguments and dis- 
cussion presented in this paper apply to the phenomenon 
of Bose-Condensation in general. 

Results from numerical simulations presented in this 
paper were carried out using typical experimental values 
for a trap containing 87 Rb atoms. The values are given 
in table |. 

TABLE I. Values of constants used in the numerical simu- 
lation. 



Constant 


Value 


m 


1.4x10"^ kg 


a 


5.2xl0" 9 m 




20.0X7T rad s" 1 
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III. SELF-CONSISTENT THERMALLY-EXCITED 
STATE GENERATION 

To overcome the prohibitive amount of computational 
time required to time evolve a thermally-excited conden- 
sate using finite temperature field theory, an alternative 
technique has been devised. This technique for generat- 
ing a self-consistent thermally-excited state involves tak- 
ing a noninteracting mean field and propagating it using 
the GP equation (Eq. (Q)), whilst the interactions are 
slowly switched on. 

One should expect the GP equation to describe the ki- 
netics and dynamics of a finite temperature boson gas 
as long as the relevant modes of excitation are suffi- 
ciently well populated. Indeed, Moore and Turok |2(| 
have shown that one can use the evolution of the classi- 
cal field during the electroweak phase transition to deter- 
mine the dynamics of the infrared modes of the quantum 
field. The classical field description neglects quantum 
fluctuations (and other quantum manifestations, such as 
squeezed states or spontaneous emission). The method 
works if, and only if, the population of the normal modes 
of the system are much greater than unity. The normal 
modes of the system can then be treated as classical ob- 
jects or quasiparticles that are statistically independent 
of the energy levels, allowing them to be represented by a 
classical field. Morgan et. al. |l^] have also stressed that 
the use of a classical field to describe a finite temperature 
BEC is valid in the limit of large mode occupations. 

The GP equation includes the dynamics and kinetics 
of a boson gas through the nonlinear term in the Hamil- 
tonian, which in effect represents all the stimulated two- 
body scattering processes. These processes ensure that 
an initial distribution will, given sufficient time, result in 
a thermal distribution of the quasiparticles, irrespective 
of the quantitative nature of the initial population dis- 
tribution. The resulting population distribution is close 
to the full Bose-Einstein distribution in the equipartition 
sense, i.e. 
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=E nm /k B T _ l 
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(10) 



where k B is the Boltzmann constant, T is the tempera- 
ture and E nm is the mode energy. In this sense the quasi- 
particles act as classical objects, macroscopically popu- 
lating the normal modes of the system. 

If the GP equation provides an accurate description of 
the quasiparticle interactions, and therefore the time de- 
pendence of the boson gas, it must also describe its tran- 
sition to equilibrium. For example, a boson gas could 
be prepared with a non-thermal distribution of the en- 
ergy level populations. Given enough time however, the 
collisional interactions within the gas would rethermalise 
it, recreating a Bose-Einstein population distribution, in 
the equipartition sense. Kagan et. al. [^0| argue that, for 
large quasiparticle populations, this is true even when 
the initial distribution is far from equilibrium or where 



the ground state is initially unpopulated. This issue has 
however, been the subject of considerable debate f23|| . 

To ensure that the resulting ensemble is fully repre- 
sentative and that the subsequent 'nonlinear mixing' is 
unbiased, the initial distribution must include all symme- 
try combinations. This requirement can be deduced from 
an explicit consideration of the parity invariance of the 
evolution of the SHO mode amplitudes. The consequent 
symmetry dependent nonlinear mixing implies that if one 
symmetry combination is omitted the resulting nonlinear 
mixing will be unbalanced. In particular, if the ensemble 
is comprised of entirely even or odd modes, resulting in 
a totally even nonlinear operator, the modes only couple 
between like symmetries producing a purely even or odd 
distribution. 

The thermalisation or nonlinear mixing of the mode 
populations described by the nonlinear term in the 
Hamiltonian of th GP equation is the underlining pro- 
cess of this method of self-consistent thermally-excited 
state generation. It allows a direct solution of the finite 
temperature analysis to be circumvented when the quasi- 
particle populations are much greater than unity. The 
self-consistent thermally-excited state generation is com- 
putationally achieved by taking a noninteracting, non- 
equilibrium, thermally-excited distribution and adiabati- 
cally increasing the nonlinear term. Adiabatically switch- 
ing on the interactions results in the relaxation of the 
mean field into an interacting, equilibrium distribution. 

To this end the initial state is a thermally distributed 
superposition of SHO modes and therefore, a thermally- 
excited solution to the N = GP equation. The initial 
state, v I'i(£,t), is given by 



o,p 



(11) 



where 0,P are the mode-cut-off numbers, \a nrn \ is the 
modulus of the mode amplitude and a nm is the mode 
phase. The SHO modes, ip n m{£,'r), in Cartesian coordi- 
nates are given by 

Vw(S,T)=if„(^)# m y^)e^ (12) 

.g-i("+ m + 1 )' r 

where H n (g) are the Hermite polynomials. 

Since the initial superposition has a finite tempera- 
ture, the initial mode probabilities are given by the lin- 
ear, /x = 0, Bose-Einstein distribution namely 
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e (n+m+l)//3 _ I 



(13) 



where the scaled temperature [3 = ksT/tuvt- The non- 
interacting thermal state has no long range order, this 
incoherence enables the mode phases of the initial state, 
&nm, to be described by a random variable on the inter- 
val (0, 27r). The numerical simulation must then be run 
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several times using different initial states to ensure that 
the kinetics are independent of the initial phases. 

The SHO modes in the initial superposition are ther- 
mally distributed, however because of the mode cut-off 
numbers, O, P, the initial mean field is a truncation of a 
thermally-excited state and therefore, not in equilibrium. 
The resulting truncation error, however, is negligible as 
the choice of initial parameters, f3 = 20 and 0,P = 20, 
ensures that the initial distribution is sufficiently accu- 
rate, whilst still being computationally manageable. Fur- 
thermore the nonlinear counterparts of the unpopulated 
SHO states, ip nm n,m> 20, should not be unfavourably 
populated as all symmetry combinations initially exist. 
Therefore the thermalisation process will couple all the 
states populating them irrespectively of the initial mode 
amplitude. 



Pre-mixing mean field density w.rt transverse coordinates 



state, Ntf, where 




Transverse distance x in SHO units 



Transverse distance y in SHO units 



FIG. 1. Surface plot of the mean field density of a typical 
initial state, with respect to the scaled spatial coordinates 
(£, 0) where O < P + 20, /3 = 20. 



The initial state, &i(£,t), is then propagated using 
the method described in section [y, where the interaction 
is adiabatically increased to its final value over the first 
half of the mixing time and the mean field is allowed to 
equilibrate over the following half. The increase in the 
nonlinearity from zero, where the initial condition is an 
exact solution to the GP equation, to its final value cor- 
responds to filling the trap with like particles, increasing 
the number of atoms in the thermally-excited state. 

To obtain a link between the number of atoms in the 
thermally-excited state and its zero temperature coun- 
terpart the final value of rj was chosen such that N was 
equivalent to the number of atoms in a Thomas-Fermi 



N TF = / |V*F(€,r)fd£ = 



4a \muJt 



1/2 



(14) 



Choosing /i = 515 results in a total number of atoms 
N = 3.07 x 10 7 . Combining this number of atoms with a 
trap ratio e = 4 x 10 9 ensures that the nonlinear factor, rj, 
is large enough to result in a nonlinear mixing rate that 
allows full thermalisation of the Bose-condensed gas. The 
large numbers involved simply reflect the slow pace of the 
nonlinear mixing in two dimensions due to the reduced 
mode density. 

The numerical method results in the thermally-excited 
state, *t(£,t to ), given by 



*T(£,r„ 



n« 
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*i(€,o) 



(15) 



where Q is the number of computational steps in the time 
dimension, St is the scaled step size and the scaled mix- 
ing time r m = QSt. The simulation was run with r m = 
100, where the nonlinearity is increased on the interval 
T m = [0, 20] and then remains constant allowing the mean 
field to equilibrate. The spatial grid has 256x256 com- 
putational points with a scaled area of 50x50 and the 
number of computational steps Q — 20000. This results 
in a calculation that requires some hours to run on a Sun 
Ultra II workstation. 

The results of the numerical simulation of the 
thermally-excited state generation are shown in Figs. |l| 
and H which show the density, p(£,r) = I'S'^, t)| 2 , of 
the thermally-excited state before and after the nonlinear 
mixing process has occurred at z = 0. The figures show 
that the thermalisation leads to a macroscopically popu- 
lated ground state, in the centre of the figure, surrounded 
by thermally-excited modes. The geometry of the initial 
state, which could have been circular but was chosen to 
be square for didactic and diagnostic purposes, has been 
changed into a circular one reflecting the symmetry of 
the trapping potential. This geometric relaxation of the 
mean field, into an equilibrium state, provides a visually 
convincing display of the nonlinear mixing. The geomet- 
ric relaxation also results in a lower temperature, a conse- 
quence of the two body scattering process acting to give 
the mean field the most probable distribution resulting 
in a maximal entropy state. 
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Transverse distance x in SHO units "' u - c Transverse distance y in SHO units 

FIG. 2. Surface plot of the mean field density of a typical 
thermally-excited state, , with respect to the scaled spatial 
coordinates (£, 0) after a scaled mixing time r m = 80. 



Further evidence of the thermalisation is shown in Fig. 
[II, which shows the populations of the normal modes for 
the initial linear and final nonlinear states for two initial 
choices of the phase, calculated using the method of adi- 
abatic following. The population distributions demon- 
strate that the nonlinear mixing is sufficient to redis- 
tribute the energy of the system into an interacting Bose- 
Einstein distribution, in the equipartition sense. This re- 
distribution was found to occur if the initial population 
distribution was Boltzmann or Bosc-Einstcin, although 
the initial choice did effect the temperature of the result- 
ing thermally-excited state. 




FIG. 3. Mode populations, \a nm \ 2 , of the initial (dashed) 
and self-consistent thermally-excited states (solid) for two ini- 
tial choices of the phase. It is plotted with respect to the 
scaled linear mode energy, n + m + 1, after a scaled mixing 
time r m = 80. 



Figs. [j-[II are typical of the results obtained, where 
the qualitative features of the final density, r m ), and 
mode populations, \a nm {T m )\ described above are inde- 
pendent of the choice of initial phases. This is shown in 
Fig. Ill and in subsequent simulations, even if the phase 
of the initial modes are chosen such that the initial state 
is coherent. 



IV. EVAPORATIVE COOLING 

With the ability to create thermally-excited states, a 
model of the evaporative cooling process can now be in- 
vestigated. The process of evaporative cooling involves 
the selective removal of the most energetic atoms in the 
thermally-excited state followed by a period of rether- 
malisation resulting in a net cooling effect. In current 
experiments the cooling is carried out by inducing a ra- 
dio frequency transition to an untrapped magnetic sub- 
level. The transition frequency is a function of the mag- 
netic field and this dependence results in the selective 
removal from the trap of atoms with a potential energy 
E > huj r f\mf\ where w r f is the angular frequency of 
the transition inducing laser. The frequency of the ra- 
diation is then gradually lowered bringing less energetic 
atoms into resonance with the untrapped state with the 
thermally-excited state rethermalising as the evaporative 
process continues. 

This energy dependent removal of the atoms makes it 
possible to view the evaporative cooling either as an ab- 
sorption of the edges of the condensate, or the lowering of 
the trap walls, allowing the escape of the more energetic 
atoms. Of the two options, it is computationally sim- 
pler to introduce an absorbing boundary, the position of 
which is a function of time. Such a model does not make 
any assumption or approximation other than represent- 
ing the state by a classical field or that the cooling can be 
modelled by a spatially dependent removal of the field. 

Inclusion of the absorption term, A(£,t), in the 
GP equation results in the time dependent nonlinear 
Schroedinger equation 



idtVfor) = ( -V 2 

^A(&r)*(£,r). 
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#(€,t) (16) 



Since the absorbing boundary must have a smooth pro- 
file to preclude the introduction of computational errors, 
A(£,t) was modelled using a Super-Gaussian 
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A(t,r) = A 1 



(17) 
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Here £ e (r) defines the position of the step function and 
Ao determines the effectiveness of the cooling which is 
related to the magnitude of the field of the evaporative 
cooling laser. The super-Gaussian defines the spatial re- 
gion that undergoes cooling and ensures a smooth tran- 
sition between the absorbing and non-absorbing regions. 
The integer m which determines the spatial width, £ s , 
over which the function steps from to Aq is set at 
20 throughout this paper, giving an estimated width of 
£ s = 0.2 x £ e (r). This choice allows Eq. (0) to rapidly 
step from to A Q whilst still varying over a sufficiently 
large range to sample on a computational grid. 

The evaporation length, £ e (r), implicitly determines 
the energy of the hotter atoms that are being removed 
from the thermally-excited state and as it is a function of 
time it also determines the dynamics of the cooling and 
ultimately the resultant ground state. 

The cooling rate depends on the collision rate and 
therefore the density of atoms in the thermally-excited 
state and it is most effective, in terms of the resultant 
condensate population, if the absorbing boundary moves 
slowly removing only the minimum number of atoms. To 
this end, two cooling strategies are compared, namely 



slightly warmer condensate than the exponential depen- 
dence. Subsequent figures presented in this paper relate 
to this linear case. 



A. Cooling: ~ size of the thermally-excited state. 



In this case the evaporation length varies monotoni- 
cally from £ e = 17 to £ e (r = 60) = 8 during the cooling 
cycle followed by 20 scaled time units of equilibration re- 
sulting in t c = 80. The change is slow enough to ensure 
that only the most energetic atoms are removed from 
the thermally-excited state. The resulting rethermalisa- 
tion process redistributes the population of the normal 
modes of the G-P equation, reducing the temperature. 
Typically the cooling cycle starts with a density, p(C), as 
shown in Fig. Ill in the previous section and results in 
the density plot shown in Fig. [I| The evidence that the 
thermally-excited state has cooled is three fold: 



(18a) 
(18b) 



where ^ is the initial evaporation length; £y is the fi- 
nal evaporation length; r c is the time scale over which 
the cooling takes place and r s = — r c / In (1 — £/ /&). The 
evap oration length has a simple linear va riat ion in Eq. 
(18a) and an exponential variation in Eq. ( |l8h| ) although 
both have identical initial and final lengths. 

Numerical simulations of the evaporative cooling pro- 
cess modelled by the a bsorbing boundary, with £ e (r) de- 
fined by Eq. (18a -181), were carried out using the nu- 
merical method described in section O. The cooled state, 
*c(£,t c ), is given by 



*c(£,t c )= J] 
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*t(£,0) 



(19) 



where the cooling time, r c = QSt. The number of com- 
putational steps and transverse grid size is the same as in 
with the cooling time, r c 



section 



III 



80, which results in 
a computation taking the same real time to complete as 
before. The initial condition used is a thermally-excited 
state, 5't(^, 0), created in the same manner as section 
III with the numerical constants remaining unchanged. 

Striking differences in the results of the evaporative 
cooling occur depending on the choice of two spe- 
cific cases are considered below. The results also show, 
however, that there is little variation between the expo- 
nential and linear dependence of the evaporation length, 
Eqs. (Il8|). Although the linear dependence, Eq. (18a), 




results in a slightly larger ground state population and a 



Transverse distance x in SH0 units - i0 Transverse distance y in SH0 units 

FIG. 4. Surface plot of the mean field density of an evap- 
oratively cooled state, $c, with respect to the scaled spatial 
coordinates (£, 0) after a scaled cooling time r c = 60 where 
£i > size of the thermally-excited state. 



a. The density of the condensate is now much 
smoother and thus the velocity, u(C) = 2mp(c) ( v ^ / *^ 1 ^ — 
vj/V^*), of the initially thermally-excited state has been 
reduced. This reduction reflects the fact that the 
thermally-excited state now has a lower temperature. 
Also the shape of the resulting condensate density has 
a profile with a greater similarity to that of the ground 
state of the GP equation. 
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Transverse distance y in SHO units 

FIG. 5. Contour plot of the phase of the mean field of an 
evaporatively cooled state, ^c, with respect to the scaled 
spatial coordinates £ after a scaled cooling time r c = 60. 



c. Fig. |6| shows the populations of the normal modes 
of the thermally-excited state before and after cooling 
for two initial choices of the phase, calculated using the 
method of adiabatic following. This figure clearly demon- 
strates that the thermally-excited state has been cooled. 
The figure also shows that the atoms from the more en- 
ergetic states have not just been removed, but that the 
population has been a transfered to the ground state as 
a result of the rethermalisation process. 



B. Cooling: £j < size of the thermally-excited state. 

In this case the evaporation length varies monoton- 
ically from £ e = 12 to £ e (Y = 60) = 7 during the 
cooling cycle, giving an initially larger spatial cut into 
the thermally-excited state. This is then followed by 20 
scaled time units of equilibration resulting in r c = 80. 
The initially larger spatial cut results in the removal of 
some atoms which may have a significant amount of ki- 
netic energy. Cooling in this fashion typically results in 
a mean field density, shown in Fig. [?], which includes a 
trapped vortex state. 
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FIG. 6. Mode populations, |a nm | 2 , of the thermally-excited 
(solid) and cooled (dash-dot) states for two initial choices of 
the phase. It is plotted with respect to the scaled linear mode 
energy, n + m+1, after a scaled cooling time r c = 60. 
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FIG. 7. A contour plot of the mean field density of an evap- 
oratively cooled state, *&c, with respect to the scaled spatial 
coordinates (£, 0) after a scaled cooling time r c = 60 where 
£i < size of the thermally-excited state. The trapped vortex 
state is on the lower left-hand edge of the condensate. 



b. The phase of the thermally-excited state, Fig. 
^, is now much smoother and has large coherent areas. 
This arises because the thermally-excited state now pre- 
dominantly consists of a ground state with only a small 
number of modestly populated excited states producing 
the phase fluctuations. 



Initially the thermally-excited state has no angular mo- 
mentum, so when viewed in an angular momentum rep- 
resentation the vortex states exist in an equal number of 
right and left handed pairs. Therefore, when the cool- 
ing is carried out as described in section 
pairs annihilate one another exactly. If, 



IV A the vortex 



however, £j is 
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less than the size of the thermally-excited state, the ini- 
tially deeper spatial cut into the thermally-excited state 
removes a vortex during the initial stages of the cool- 
ing cycle resulting in an unbalanced vortex population. 
Consequently, as the thermally-excited state is cooled, 
at least one vortex will remain trapped with the ground 
state, resulting in a state with a non zero angular mo- 
mentum. 

The resulting trapped vortex state is robust and re- 
mains with the condensate until the spatial confinement, 
due to the evaporative cooling, causes the state to be 
energetically unfavourable at which point it disappears. 

This method of vortex trapping may be directly appli- 
cable to two-dimensional experimental situations if the 
evaporation surface in the experiment can be setup so 
that it removes atoms with a proportion of kinetic en- 
ergy in the initial stages of cooling. 



V. CONCLUSION 

We have shown in this paper that using the GP 
equation to time evolve an initial state is a numeri- 
cally efficient method for the creation of self consistent 
thermally-excited states. Results presented demonstrate 
that the GP equation acts to thermalise an initially non- 
equilibrium state through two body scattering processes 
described by the nonlinearity within it. The resulting 
nonlinear mixing generates a self-consistent thermally- 
excited state with an interacting Bosc-Einstcin popula- 
tion distribution, in the equipartition sense. Where the 
ensemble average of the GP equation describes the finite 
temperature kinetics and dynamics of a boson gas if and 
only if the quasiparticle populations |a nm | 2 ^> 1. 

A model of the evaporative cooling process without the 
assumption of sufficient ergodicity or approximations re- 
lating to the strength of the interactions between the 
thermal and condensate parts of the thermally-excited 
state was presented. This model uses an extension of the 
GP equation based on the removal of the most energetic 
atoms in the spatial distribution. Results from a two- 
dimensional numerical simulation show that the evapo- 
rative cooling model produces a visibly cooler mean field 
density and a population transfer to the ground state, 
which gives rise to a coherent condensate. These results 
demonstrate that the model accurately reflects the evap- 
orative cooling process and consequent rethermalisation 
of the thermally-excited state. Furthermore other results 
presented show that under a different cooling regime vor- 
tex states can be trapped along with the condensate and 
it is noted that it may be possible to achieve this using 
current experimental configurations. 

The possibility of creating self-consistent thermally- 
excited states and the numerical simulation of these 
states in two dimensions allows the study of many 
new and interesting phenomenon such as vortex and 
other quasiparticle interactions, coherence properties of 



thermally-excited states and more accurate models of 
Bose-condensation experiments. 
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